# This file produces tables 7 and 8 and figure 2
set.seed(123)
# Set working directory
setwd("/Users/ericguntermann/Documents/Papers/Representation of Party Preferences/Replication")

## Criterion 3 (proportional)
criteria <- read.dta("criteria.dta")
criteria <- criteria[-which(criteria$countryyear %in% c("Thailand 2001", "Thailand 2011","Japan 1996", "Japan 2013")),]
d1 <- with(criteria, list(lddiff=lddiff, proportional=proportional, gdppercap=gdppercap, freedomhouse=freedomhouse))
d1$N <- length(unique(criteria$countryyear))

m1 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.proportional*proportional[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.proportional ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

p1 <- c("b", "b.proportional", "b.gdppercap", "b.freedomhouse")

i1 <- function(){
  list("b"=rnorm(1), "b.proportional"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}



## Criterion 3 (gallagher)
library(foreign)
dat <- read.dta("criteria.dta")
d2 <- with(dat, list(lddiff=lddiff, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse))
d2$N <- length(unique(dat$countryyear))
m2 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i2 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}
p2 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse")


## Criterion 3 (mdm)
library(foreign)
criteria <- read.dta("criteria.dta")
d3 <- with(criteria, list(lddiff=lddiff, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse))
d3$N <- length(unique(criteria$countryyear))

m3 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p3 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse")


i3 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}

## Criterion 3 (mdm) with number of parties
library(foreign)
criteria <- read.dta("criteria.dta")
d4 <- with(criteria, list(lddiff=lddiff, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d4$N <- length(unique(criteria$countryyear))
m4 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i] + b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p4 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse", "b.ngov")


i4 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}


d <- list(d1,d2,d3,d4)
p <- list(p1,p2,p3,p4)
m <- list(m1,m2,m3,m4)
ini <- list(i1,i2,i3,i4)


library(foreach)
library(doParallel)
library(R2jags)
registerDoParallel(4) 
cl <- makeCluster(4)
registerDoParallel(cl)
results <- foreach(i=1:4, .packages = c('R2jags')) %dopar% jags(data=d[[i]], model=m[[i]], parameters.to.save=p[[i]], inits=ini[[i]], n.iter=100000, n.burnin=50000, n.thin=10, n.chains=2)


results1 <- as.mcmc(results[[1]])
results2 <- as.mcmc(results[[2]])
results3 <- as.mcmc(results[[3]])
results4 <- as.mcmc(results[[4]])

# Convergence diagnostics

library(ggmcmc)
p1 <- ggs_geweke(ggs(results1)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 1")
p2 <- ggs_geweke(ggs(results2)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 2")
p3 <- ggs_geweke(ggs(results3)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 3")
p4 <- ggs_geweke(ggs(results4)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 4")
library(gridExtra)

pdf("criterion3_conv.pdf")
grid.arrange(p1,p2,p3,p4)
dev.off()

## GET OUTPUT
#Model 1 (Proportional)
coefs1 <- as.data.frame(as.matrix(results1))
b <- coefs1[,"b"]
b.proportional <- coefs1[,"b.proportional"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co1 <- c(mean(b), mean(b.proportional), NA, NA, mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.proportional), NA, NA, sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.proportional)>0, mean(b.proportional>0),mean(b.proportional<0)), NA, NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 2 (Gallagher)
coefs1 <- as.data.frame(as.matrix(results2))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]


co2 <- c(mean(b), NA, mean(b.gallagher), NA, mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), NA, sd(b.gallagher), NA, sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 3 (MDM)
coefs1 <- as.data.frame(as.matrix(results3))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co3 <- c(mean(b), NA, NA, mean(b.mdm), mean(b.freedomhouse), mean(b.gdppercap))
sd3 <- c(sd(b), NA, NA, sd(b.mdm), sd(b.freedomhouse), sd(b.gdppercap))
p3 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))

# Table 7
library(xtable)
model.effects <- data.frame(co1,sd1,p1,co2,sd2,p2,co3,sd3,p3)                                                                                                                
rownames(model.effects) <- c("Intercept", "Proportional", "Gallagher*", "MDM*", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P","Mean", "SD","P")
xtab <- xtable(model.effects, caption="Criterion 3", digits=c(2))
sink("criterion3.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)


#Model 4 (MDM with ngov)
coefs1 <- as.data.frame(as.matrix(results4))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]


co1 <- c(mean(b), mean(b.mdm), mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b), sd(b.mdm), sd(b.ngov), sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)),ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


# Table 8
library(xtable)
model.effects <- data.frame(co1,sd1,p1)                                                                                                                
rownames(model.effects) <- c("Intercept", "MDM*", "Number of parties", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 3", digits=c(2))
sink("criterion3b.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)


# Make Figure 2

library(ggplot2)


library(foreign)
library(R2jags)
set.seed(123)
criteria <- read.dta("criteria.dta")
d2 <- with(criteria, list(lddiff=lddiff, gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d2$N <- length(unique(criteria$countryyear)) + 600
d2$ngov <- c(d2$ngov, rep(seq(1,6),100))
d2$gdppercap <- c(d2$gdppercap, rep(mean(d2$gdppercap),600))
d2$freedomhouse <- c(d2$freedomhouse, rep(mean(d2$freedomhouse),600))
d2$lddiff <- c(d2$lddiff, rep(NA, 600))
m2 <- function(){
  for(i in 1:N){
    lddiff[i] ~ dnorm(mu[i], tau)
    mu[i] <- b + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i] + b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p2 <- c("b","b.gdppercap", "b.freedomhouse", "b.ngov", "mu")

i2 <- function(){
  list("b"=rnorm(1),"b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}
results.pred <- jags(data=d2, model=m2, parameters.to.save=p2, inits=i2, n.iter=100000, n.thin=100, n.burnin=20000, n.chains=2)
results.mcmc <- as.mcmc(results.pred)
coefs1 <- as.data.frame(as.matrix(results.mcmc))
b <- coefs1[,"b"]

b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]
library(gtools)
mu <- coefs1[,grep("mu", colnames(coefs1))]
mu <- mu[,mixedsort(colnames(mu))][,88:687]

df <- data.frame(ngov=criteria$ngov, y=criteria$lddiff, PR=criteria$proportional)
df$cy <- criteria$countryyear
df$PR[df$cy %in% c("Japan 1996", "Thailand 2001", "Thailand 2011")] <- 3
df$PR[df$PR==1] <- 2
df$PR[df$PR==3] <- 1
df$PR <- as.factor(df$PR)
levels(df$PR) <- c("Non-PR","Mixed Parallel", "PR")
df <- df[order(df$y),]
df.id <- df[c(1:3,85:87),]
colnames(df.id) <- c("xngov", "yld", "pr", "cy2")

df2 <- apply(mu, 2, quantile, probs=c(0.05, 0.5,0.95))
df2 <- data.frame(ngov=1:6, t(df2))
colnames(df2) <- c("ngov", "lwr", "y", "upr")

library(ggplot2)
p <- ggplot(df, aes(x=ngov, y=y)) + geom_point(aes(shape=PR, colour=PR), size=2) + geom_ribbon(data=df2,aes(ymin=lwr,ymax=upr), alpha=0.3) + geom_line(data=df2,aes(x=ngov, y=y)) + geom_hline(data=df, aes(yintercept=mean(y)), alpha=0.9, linetype="dashed") + xlab("Number of parties in government") + scale_colour_grey(start = .6, end = 0) + theme_bw() + ggtitle("") + ylab("Criterion 3 scores") + theme_set(theme_bw(base_size = 14)) + geom_text(data=df.id,aes(x=xngov, y=yld + 0.1, label=cy2), size=4) + geom_text(aes(0.5,mean(df$y),label = "Mean=0.99"), vjust=-0.5, hjust=0.75, colour="black", size=5) + xlim(-.5,6.5) + theme_bw()
ggsave("criterion3plotngov.pdf")
